Circuit QED with fluxonium qubits: theory of the dispersive regime 
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In circuit QED, protocols for quantum gates and readout of superconducting qubits often rely on the dispersive 
regime, reached when the qubit-photon detuning A is large compared to the mutual coupling strength. For 
qubits including the Cooper-pair box and transmon, selection rules dramatically restrict the contributions to 
dispersive level shifts \. By contrast, in the absence of selection rules many virtual transitions contribute to \ 
and can produce sizable dispersive shifts even at large detuning. We present theory for a generic multi-level 
qudit capacitively coupled to one or multiple harmonic modes, and give general expressions for the effective 
Hamiltonian in second and fourth order perturbation theory. Applying our results to the fluxonium system, 
we show that the absence of strong selection rules explains the surprisingly large dispersive shifts observed in 
experiments and leads to the prediction of a two-photon vacuum Rabi splitting. Quantitative predictions from 
our theory are in good agreement with experimental data over a wide range of magnetic flux and reveal that 
fourth-order resonances are important for the phase modulation observed in fluxonium spectroscopy. 

PACS numbers: 85.25.Cp, 03.67.Lx, 42.50.Pq 



I. INTRODUCTION 

In many respects, the quantum physics of superconducting 
circuits 1-3 resembles that of atoms: both systems feature a 
set of discrete, non-equidistant energy levels, and both can 
be probed by virtue of their interaction with photons. Within 
circuit QED 4-6 , this interaction is harnessed to manipulate and 
measure the quantum state of the superconducting circuit with 
great success. 

In many realizations of the circuit QED architecture, the 
dispersive regime plays a particularly important role for im- 
plementing the readout and gate operations required for a uni- 
versal quantum computer. The general idea behind the disper- 
sive regime is simple: when detuning the qubit frequency far 
from the frequency of the resonator, the interaction-induced 
conversion of a qubit excitation into a photon becomes inef- 
fective. Specifically, the probability amplitude for the con- 
version is proportional to the small parameter gj A, where g 
denotes the coupling strength and A the detuning between 
the qubit and photon frequency. Accordingly, the dispersive 
regime is the primary setting for performing qubit gates. 4 ' 7 ~ n 

In the dispersive regime, the qubit-photon coupling man- 
ifests in the form of energy shifts of Lamb and ac-Stark 
type. 4 ' 7 ' 8 ' 12-19 Here, ac-Stark shifts correspond to state- 
dependent energy shifts which, in the simplest case, take on 
the form %a i ^ aa z . This expression can be interpreted as a fre- 
quency shift of the resonator, with the size of the shift depend- 
ing on the state of the qubit, or alternatively, as a shift of the 
qubit transition frequency, with the size of the shift depend- 
ing on the photon state of the resonator. Consequently, the 
dispersive regime provides both a convenient means of qubit 
readout. 7 ~ 9 ' 12 ~ 15 ' 19 ~ 22 , as well as a measurement tool for in- 
vestigating the resonator state. 12,23 Possible limitations to this 
simple picture due to corrections of higher order in the pa- 
rameter g/A have recently been studied by Boissonneault et 

al 16,17 

The physics of the dispersive regime becomes richer when 



higher levels of the superconducting circuit (which we hence 
refer to as qudit) participate in the virtual transitions that con- 
tribute to the dispersive shifts. The simplest manifestation of 
this is the contribution of the third level of the transmon to the 
dispersive shift, even under conditions when real occupation 
of this level is negligible. For specific level configurations 
relative to the resonator frequency, the two partial contribu- 
tions to x a dd up constructively and give rise to the straddling 
regime with characteristically large dispersive shifts. 24 - 25 Re- 
cently, enhanced dispersive shifts have also been predicted 
and confirmed experimentally for a flux qubit coupled to a 
resonator. 26 Similar to the transmon straddling regime, higher 
levels of the flux qubit are responsible for the observed shift 
enhancement. 

In this paper, we present theory that systematically de- 
scribes the dispersive regime for a generic circuit QED sys- 
tem consisting of a multi-level qudit coupled capacitively to 
one or multiple harmonic modes. In Section II, we derive the 
general expression of the effective Hamiltonian governing the 
dispersive regime up to (and including) terms of fourth order 
in the coupling between the qudit and the harmonic modes. 
We verify that our results correctly reproduce the well-known 
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FIG. 1 : (Color online) (a) General scheme of circuit quantum electro- 
dynamics with capacitive coupling between one (or several) resonant 
mode(s) and a superconducting circuit acting as a qudit. Examples 
of superconducting qudits are atoms are (b) the Cooper pair box in 
the charging and transmon regime, and (c) the fluxonium circuit. 



2 



expressions for the dispersive regime of the Cooper pair box 
in the charging and transmon regime. In Section III we then 
apply our results to the fluxonium system 27-29 where, different 
from the Cooper pair box case, the lack of selection rules al- 
lows for a large number of terms to contribute to the dispersive 
shifts. We compare our theoretical predictions with the data 
from reflection and spectroscopy experiments obtained previ- 
ously by the Yale group, 27 ' 29 ' 30 and summarize our findings 
and conclusions in Section IV. 



cjj-.w, Equations (1) and (2) may model, for instance, circuit 
QED systems based on the Cooper pair box, the transmon, 24 ' 35 
or the recently developed fluxonium device. 27 All three ex- 
amples are summarized in Table I. Note that the interaction, 
in general, does not conserve the overall excitation number 
Si a l a j + 'I ( ' I un l ess special selection rules restrict 
the coupling cjj-iv to nearest-neighbor qudit levels and rotat- 
ing wave approximation (RWA) is assumed. 



II. DISPERSIVE REGIME OF CHARGE-COUPLED 
CIRCUIT QED SYSTEMS 

A. General Model 

We consider a charge-coupled circuit QED system 4-6 as 
depicted in Fig. 1. The system contains a superconducting 
circuit 1 ' 3 ' 31 acting as a quantum system with a discrete and 
anharmonic energy spectrum. Honoring the multi-level na- 
ture, we will refer to it as a qudit and denote its eigenstates 
and corresponding energies as 1 1) and e;, respectively. Here, 
/ enumerates the eigenstates starting with I = for the qu- 
dit ground state. The qudit is coupled to microwave photons 
inside one or several 32 superconducting transmission line res- 
onators) with normal mode frequencies ojj. We consider a 
finite number of such modes enumerated by j = 1,2,... The 
inclusion of the infinite set of higher modes inside each res- 
onator is beyond the scope of the present paper, but has been 
subject of recent studies 33 34 indicating that the qubit size pro- 
vides a natural cutoff in j. 

The generic model Hamiltonian H = Hq + V describing 
the coupled circuit QED system thus captures the bare qudit 
and harmonic modes, 

H =^2u j a]a J +^2ei\l)(l\ (1) 
j i 

and the interaction between them, 

y =EE^IO(^'l(« J +«])- (2) 

Here, aj (aj-) is the usual annihilation (creation) operator for a 
photon in mode j. The interaction term describes the coupling 
between the relevant charge variable of the qudit, (2e)IM, and 
the electric voltage of the resonator at the qudit position, Vj — 

VJ ms (a'j + aj). The resulting coupling coefficients are given 
by 

g rM , =g J (l\N\l / ), (3) 

where gj = 2e/3 ; V™ ns abbreviates the qudit-independent part 
of the coupling strength and f3j is a dimensionless capacitance 
ratio, typically of order unity. 24 We set h = 1 and, in the 
following, only deviate from this convention when discussing 
concrete experimental parameters and observables. 

The above Hamiltonian provides a generic model for sim- 
ple, charge-coupled circuit QED systems. Adopting the ap- 
propriate qudit energy spectrum {e;} and coupling parameters 



B. Effective Hamiltonian for the Dispersive Regime 

The qudit-photon interaction V facilitates transitions | / } — > 
| ) between qudit states which are accompanied by the emis- 
sion or absorption of photons. The dispersive regime of cir- 
cuit QED 4 ' 16-18 ' 36 describes the situation when such transi- 
tions are suppressed due to large detuning between the qudit 
transition frequencies and the relevant mode frequencies. The 
adiabatic elimination of the coupling V, appropriate when the 
energy mismatch is large compared to the coupling strength, 
has found many applications in different branches of physics 
and, depending on context, is known under several names in- 
cluding van-Vleck perturbation theory 37 and Schrieffer- Wolff 
transformation. 38 

Adiabatic elimination is illustrated most simply for an en- 
ergy spectrum featuring a large gap between two groups of 
unperturbed states. By construction of an appropriate canon- 
ical transformation, one obtains an effective Hamiltonian in 
which the weak interaction of states above the gap with states 
below, is eliminated in favor of dressed states with energies 
slightly shifted relative to the unperturbed ones. The setting 
of two subspaces separated by a single large gap, however, 
is not a necessary requirement for the approach. The typical 
situation of the dispersive limit in circuit QED indeed differs 
from that simple setting: each unperturbed state | ni ) o forms 
its individual subspace as long as transitions from one qudit 
state to another remain sufficiently detuned from the photon 



TABLE I: Level structure and selection rules for three different su- 
perconducting (sc) qubits, capacitively coupled to a single resonator 
mode. The frequency cu p = ^J&EjEc denotes the plasma oscil- 
lation frequency, and Ej and Ec are the Josephson and charging 
energy, respectively. For the Cooper pair box in both the charging 
(Ej <5cC Ec) and transmon (Ej 3> Ec) regime, simplifying selec- 
tion rules apply. No simple selection rules exist for the fluxonium 
device. (Abbreviations used: TLS - two-level system, MLS - multi- 
level system.) 
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frequencies. Here, stats are labeled by the set of photon num- 
bers n = (ni,n2, • ■ •) and qudit states I = 0, 1, Idots The 
contribution of several virtual transitions to each energy level 
shift can make the physics of the dispersive regime quite rich. 
In part, this is already true for transmon-based circuit QED 
systems, where the I = 2 state gives rise to enhanced level 
shifts in the straddling regime. 24 Even more interesting fea- 
tures emerge in the dispersive regime of the fluxonium de- 
vice, which we discuss in detail in Section III. Based on the 
lucid description given by Cohen-Tannoudji et al., 39 we next 
summarize the systematic procedure for obtaining the effec- 
tive dispersive Hamiltonian of the circuit QED model speci- 
fied in Eqs. (1) and (2). 

As a necessary condition for the validity of the dispersive 
approximation, all one-photon transitions among low-lying 
qudit states must be strongly detuned from the harmonic mode 
frequencies. To formulate this condition quantitatively, we in- 
troduce compact notation for transition energies and detun- 
ings: t\y = e i — e i' abbreviates the energy released in the 
qudit transition I— >l' (note that ew is negative when I' > I), 
and Aj;K' = ew ~ u>j denotes the detuning between this tran- 
sition and resonator mode j. In this notation, the condition for 
the dispersive regime reads 

l A i;M'l » \9i;U'W n i + 1 , ( 4 ) 

where the photon number is typically restricted to rij = 
when assuming dilution refrigerator temperatures fcgT -C cjj, 
but may reach higher values when the system is driven with 
microwave tones. 

Condition (4) motivates the perturbative treatment of the 
interaction V which couples the unperturbed Hq eigenstates 
| nl )o. 46 Each eigenstate of H is a dressed state with the ma- 
jority of all probability amplitudes in a single state | nl ) . As 
a result, the labeling of bare states can be maintained for the 
dressed eigenstates | nl ) = e~~ lS \ nl )q. The diagonalization 
of H, up to a specified order in V, is achieved by the unitary 
transformation H' = e He~ lS . Note that, by construction, 
eigenstates of H' are just the unperturbed states | nl )q. 

The procedure 39 for obtaining the required hermitean gen- 
erator S now follows from the two conditions that H' be di- 
agonal and the generator S be off-diagonal in the unperturbed 
basis: 

PnlH'P^y ~6 nn ,6 w and P n iSP nl =0. (5) 

Here, P n i = \ nl )oo( n£ | is the projector onto a single unper- 
turbed state. Introducing an auxiliary parameter A for count- 
ing powers in V, one constructs S — J2m=i ^ m Sm and 
H' = H + X)m=i ^ m -^m order by order, by comparing with 
the nested commutator series 

oo 

H' = e lS He~ lS = V — [iS, H] m (6) 
^— ' m! 

and enforcing the conditions (5). (For further details, see Ap- 
pendix A.) 

In practice, this procedure quickly becomes cumbersome 
for terms beyond second order. Since terms of fourth order in 



the interaction will turn out to be relevant for the dispersive 
regime of fluxonium, we devise a way to bypass the evalu- 
ation of fourth-order nested commutators as follows. Note 
that there is a natural equivalence between the construction of 
the generator S on one hand, and the ordinary form of time- 
independent perturbation theory (yielding corrections to ener- 
gies and states) on the other hand. The two approaches merely 
differ in whether the basis change to the approximate eigen- 
basis is carried out as an active or passive transformation. In 
the first approach, the perturbation determines the generator 
S which brings the Hamiltonian H into the diagonal form 
of H' . In the second approach, the perturbation affects the 
dressed states which one constructs explicitly in the unper- 
turbed basis as \ nl) = e~ lS \nl)o. Fortunately, obtaining 
higher-order corrections for eigenenergies in the second ap- 
proach generally does not involve nested commutators. Thus, 
we first establish the generic structure of H ' up to the desired 
order, leaving all energy coefficients of individual terms to be 
determined. We then apply the inverse unitary transformation 
(cut off at the same order) and obtain the effective Hamilto- 
nian H e ff = e~ lS H' e lS . Finally, we employ ordinary pertur- 
bation theory to find the eigenenergy corrections and extract 
from them the undetermined energy coefficients to complete 
the effective Hamiltonian. 

The generic form of H' is dictated by the conditions from 
Eq. (5), and is easily obtained as follows. Since Eq. (5) ex- 
cludes all coupling between different subspaces, H' can be 
expressed as 

H'=J2 PmH'P nl = Y, E niPni- (7) 

n,/ n,l 

Any operator contributing to H' = Y2k H'k must be diagonal 
in the unperturbed basis, i.e., P n iH' k P n i ^ 0. Evidently, each 
such contribution can only consist of harmonic-mode number 
operators and qudit projectors. The resulting general form, 
after performing the inverse unitary transformation, is 

H em =a k \{{^ ] ) N ^\l k ){l k \. (8) 

j 

Here Nj k > are integer exponents, a k is an energy coeffi- 
cient, and the harmonic oscillator and projection operators are 
dressed-state operators, i.e., 

•'J, <'"<',, ' s and \l k )(l k \ =e lS \l k ) ao (lk\e- %s . 

(9) 

Since the interaction V is of order one and consists of a sum 
over operator terms with only one harmonic ladder opera- 
tor each, the perturbation order A m of any contribution (8) 
cannot be smaller than the number of ladder operators, i.e, 
m > 2 Y^j Nj k . This provides us with the necessary informa- 
tion to obtain the generic structure of the effective Hamilto- 
nian. 

Second-order terms. — The generic structure of the effective 



4 



Hamiltonian in second-order perturbation theory is 

H eB =5^c; i ata i + 5^e,|0<i| 

+ X> i iata i |0<i|+X>IO<*l- 

3,1 I 



(10) 



n + ej,li 



n,l ' ■> 



n — e,-, (j 



Here, the third and fourth terms describe dispersive shifts of 
ac-Stark type and qudit level shifts of Lamb type. As usual, 
the ac-Stark shifts may manifest as harmonic-mode frequency 
shifts which depend on the occupied qudit level I, ojj — > 
LJj + Xj;l> or as qudit energy shifts which depend on photon 
numbers, e; —> e; + J^j Xj-,i a ] a j- 12 ' 23 Note that terms of the 
form Awj-ajaj corresponding to pure shifts of resonant mode 
frequencies may be absorbed by letting Xj-l Xj;l + Auj j- 
To determine the coefficients an d K i we use the ordinary 
expression for the second-order energy correction: 



E. 



(2) 
III 



nl\V\n'l') o0 (n'l'\V\nl] 



E, 



(0) 



E. 



(o) 



(11) 



where the primed sum indicates that the term n'l' = nl is to 
be omitted. For further evaluation, we separate the interaction 
into photon creation and annihilation terms of the individual 



modes, V = EjOj + ) wnere 



LI' 



and V; 



(12) 



The product of transition matrix elements in the numerator of 
Eq. (11) selects the combinations V^~V~ and Vj~Vj + , which 
change the photon number by one and subsequently undo this 
change. The virtual transitions affecting photon number and 
qudit levels are illustrated conveniently in the ladder diagram 
shown in Fig. 2. They yield 

E nl = Ej,J'K-[Xi;H' - Xjd'l] + Xj;ll') 
for the second-order energy correction, where 

Xy,w = \9r,w \ 2 / A 3-,w d3) 

abbreviates partial dispersive shifts. The wanted energy coef- 
ficients in Eq. (10) can now be read off. The resulting expres- 
sions are given by 



Xjd 



^2(Xj;W - Xj;l'l) 



and 



XM 



(14) 



3,1' 



Note that both expressions include a summation over all qudit 
levels Thus, higher qudit levels - even when unoccupied - 
may contribute substantially to the dispersive shifts of photon 
frequencies and lower qudit levels. This fact is also clearly 
illustrated with the level diagram shown in Fig. 3. 

We have verified that the same expressions are obtained 
with the active transformation method (see Appendix A). In 
addition, one can readily confirm that Eqs. (10) and (14) cor- 
rectly reproduce the following well-known results. First, for 



FIG. 2: (Color online) Second-order ladder diagram showing the 
contribution of harmonic mode j to the eigenenergy correction for 
state | n, I ). The two interfering paths correspond to virtual transi- 
tions which intermediately decrease (K j + V ) _ ) and increase (V~ V^ + ) 
the photon number in mode j. Without selection rules, summation 
includes all qudit levels Zi. ej denotes the unit vector with "direc- 
tion" j. 
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FIG. 3: (Color online) Level diagram illustrating the second-order 
dispersive shifts for a multi-level qudit coupled to a single photon 
mode. For each column (definite photon number), solid lines show 
the bare energy levels with I — 0, 1 and 2 qudit excitations. Diagonal 
dashed lines show the relevant couplings between the I — 0, 1 levels 
with n a photons (middle column) to states with n a ± 1 photons (red 
lines: coupling I — 0, 1 states, blue lines: coupling to higher states). 
Note that higher levels (e.g. 1 = 2) participate in the couplings and 
affect the dispersive shifts. The shifted levels for the I = 0, 1 states 
are shown as horizontal dashed lines. In this particular example, the 
qudit transition 621 is nearly resonant with the photon frequency uj a 
leading to a large shift of the I = 1 level. 



the Jaynes-Cummings model with a two-level system (TLS) 
coupled to a single resonator mode we have gi-w = g for 
(I, I') — (1, 0) or (0, 1). For all other choices, gx-u> vanishes. 
We thus obtain 



H. 



eff; JC 



eio + X01 



°* + Xoi a acr z + const. (15) 



which agrees with the known result. 4 Second, we consider the 
multi-level transition device coupled to one resonator mode. 24 
In this case, there is a selection rule allowing only for coupling 
of nearest-neighbor transmon levels, i.e. g\-iv — for V ^ 
I ± 1. With this, we recover 

#eff; transmon = + ^ e l I 1 ) ( 1 I ( 16 ) 

I 

+ l>J,*-i - Xi+i, ; )a f a| + XU-i\ ( ^1, 



2>0 
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which correctly leads to the expression (xoi — Xi2/2)a^aer 2 
for the ac-Stark term upon projection onto the subspace 
spanned by the I = and 1 qudit states. 24 

Fourth-order terms. — The fourth-order terms of the Hamil- 
tonian take the form 

E X» 



H eS = [Eq. (10)] 



E 

3,1 



3,1 
) 2 \l){l\ 



E 

I 



E ^ a l a i a J a jlO( / l 



(17) 



Beyond additional corrections to terms already present in sec- 
ond order [first line of Eq. (17)], the fourth-order terms also 
introduce interaction among harmonic modes of self-Kerr and 
cross-Kerr type. We denote the corresponding coefficients by 
rjj-i and respectively. For simplicity (and also motivated 
by the experimental data to be discussed in Section III), we fo- 
cus our discussion on a set of two harmonic modes and refer 
to them as a and b mode. 

The fourth-order corrections to the eigenenergies are given 
by 40 



E, 



(4) 



E 

M,P,Q 



' ^NM^pVpqVqn 



v ' |^nm| 2 w|Kp| 2 



(18) 



Here, N, M, P and Q are multi-indices of the form {n, I}. En- 
ergy differences in the denominators are given by E NU = 
E^ — E^\ Matrix elements are abbreviated by V NN ' = 
(n| V|n'). Terms involving diagonal matrix elements have 
been dropped in Eq. (18) since V NN = in the cases of in- 
terest. We next sketch the evaluation of the fourth-order cor- 
rections and provide the relevant ladder diagrams. 

We denote the two terms on the right-hand side of Eq. (18) 
by (i) and (II) and start with discussing the latter. Term (II) is 
the product of two factors with a structure nearly identical to 
the second-order expression (11) and thus evaluates to 



3,h 

X E^'tXj'^ -Xj';hl] +Xj':lh)- (19) 

3',h 

By expanding this expression, we identify contributions to 
each of the fourth-order terms given in Eq. (17). In addition to 
the simple poles like (e«/ —ujj)^ 1 which already appear in the 
second-order energy coefficients, new double and triple poles 
with the same denominators (ew — uij) emerge. 

Next, we turn to term (i) in Eq. (18) which cannot be factor- 
ized. We classify the contributions from (i) in terms of ladder 
diagrams. The rules governing these ladder diagrams are as 
follows: 

1. Each ladder step is labeled by a multi-index (n m , l m ) 
specifying the occupation of the harmonic modes and 
the qudit level. 

2. Starting from the right, each subsequent ladder step is 
related to the previous one by a virtual transition ef- 
fected by the operator V- (label on the arrow). 



3. The set of all possible paths is constrained by the condi- 
tion that the left-most state must coincide with the right- 
most state (n, I). Thus, each path must contain as many 
V+asVf. 

4. Each path traversing from right to left gives a contri- 
bution involving summation over all intermediate qudit 
levels h,h,h- 

5. The product of matrix elements in the numerator of each 
term is determined by the sequence of arrow labels in 
each path. For example, the sequence V^V 2 ~V 2 + V{ h 
results in the product ffi^S^WaS^WiSi;"!- 

6. The denominator of each term consists of a product of 
three energy differences of the form E^ 1 — E^X where 
(n', I') labels the virtual intermediate states on the inner 
three ladder steps. 

The extraction of the various energy coefficients can fur- 
ther be simplified by noting that the ladder diagrams also 
directly specify the operator structure of the resulting terms 
in the effective Hamiltonian. For example, the numera- 
tor ViV^V^V^ is associated with terms of the structure 
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We first treat the single-mode contributions shown in Fig. 
4(a). For each given mode j, the diagram gives rise to six 
different paths resulting in terms of the form 



E E 

l ^1,^2)^3 



9j;ll39j;hh9j;hh9 j;hl , 
E3E2E1 



op(a J -,a})|Z><i|. (20) 



The corresponding energy denominators and operators are 
specified in the table accompanying Fig. 4. The first and 
fourth terms in Fig. 4(b) exhibit new poles absent in second- 
order perturbation theory. In frequency space, these poles 
occur when the conditions e« 2 ± 2u)j = are met and thus 
signal additional resonances when qudit transitions match the 
energy of two photons in mode j. We will argue below 
that such resonances have indeed been observed in previous 
experiments. 27 ' 29 ' 30 

With this motivation in mind, we turn to the remaining 
fourth-order contributions shown in Fig. 5. All of them are 
of dual-mode type, i.e., they include participation of two dif- 
ferent harmonic modes in the virtual transitions. For clarity, 
we label the two modes j = a and j' = b 7^ a in the follow- 
ing. All resulting contributions can be cast into the form 
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E 



94939291 
E3E2E1 



op(a,a+,b,bt)|0(/|. (21) 



where g v = gj v ;i v \ v _ x with j v G {a, b} and l 4 = l = I. 
The harmonic lowering operators for the two modes are now 
simply denoted by a, b. All the possible paths are listed in 
Fig. 5, and the corresponding coefficients and operators are 
summarized in Fig. 5(d). 

New types of poles arise in the ladder diagrams of Fig. 5(b) 
and (c) when one of the resonance conditions 
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FIG. 4: (Color online) (a) Single-mode ladder diagram involving 
only virtual excitations of mode j. (b) The corresponding table spec- 
ifies terms according to Eq. (20). An additional prime on E 2 signals 
exclusion of the term l 2 = I from the sum over l 2 . 



is met. Such resonances occur whenever a qudit transition 
matches either the energy required for placing one photon 
each in mode a and mode b, or the energy required to con- 
vert an a photon into a b photon. We find that these additional 
resonances lead to observable effects and can be pinpointed in 
the data from previous experiments, 27,29 ' 30 as we will discuss 
in the following section. 



III. APPLICATION: DISPERSIVE REGIME OF THE 
FLUXONIUM DEVICE 

Equipped with the general expressions for the effective 
Hamiltonian, we now study their concrete application to the 
fluxonium circuit. 27 The fluxonium device is of particular in- 
terest in this context since experiments have shown surpris- 
ingly large dispersive shifts 29 and, as we will see, transitions 
between fluxonium states are not strongly restricted by se- 
lection rules. The simplest model of fluxonium is obtained 
by shunting a small Josephson junction with the large kinetic 
inductance from a Josephson junction array [see Fig. 1(c)]. 
For fluxonium, the flux-dependent eigenenergies ei(§ ext ) and 
corresponding eigenstates {| I)} are thus determined by the 
Hamiltonian 

H t = 4£ C N 2 - Ej cos((^ - 2vr$ ext /$o) + \e l ^ (22) 

Here, the operator N = Q/2e describes the charge on the 
junction capacitance (in units of the Cooper pair charge) . Its 
conjugate operator 99 = 27r<£>/$o describes the loop flux in 
such units that ip coincides with the phase difference across 
the junction. Circuit quantization imposes the commutation 
relation [<p, N] = i between the two operators. Analogous to 
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FIG. 5: (Color online) (a)-(c) Dual-mode ladder diagrams contribut- 
ing to fourth-order eigenenergy corrections. For a given selection of 
two modes, labeled a and 6, each diagram has an associated diagram 
obtained by interchanging the roles of a and b. Only for (a) one finds 
that label interchange (a o b) yields an identical expression for the 
resulting contribution, (d) Table detailing the resulting contributions 
according to Eq. (21). 



the usual relation between position and momentum operators, 



' dip 



in if 



the charge operator N thus takes the form N = 
representation. The Hamiltonian Hf is accompanied by the 
usual boundary condition \im ip ^± OD ipi(<p) = 0. 

While analytical expressions can be obtained in the limit 
of very small inductive energies i?L, 28 numerical diagonal- 
ization remains the most useful approach for the intermediate 



parameter values realized in the experiments. 27 ' 29 ' 41 Employ- 
ing the harmonic oscillator basis (which diagonalizes Hf for 
Ej = 0) and taking the necessary precautions for conver- 
gence with respect to the basis truncation, we obtain the en- 
ergy spectrum and eigenstates and calculate the charge matrix 
elements (l\ N | Z' ) entering the coupling parameters gj-w- 
For the interested reader, we provide details of the numerical 
diagonalization scheme in Appendix B. 

Figure 6 shows representative results for the charge matrix 
elements obtained in this way, using model parameters which 
match the experimental values. 41 For the example of magnetic 
flux $ext = 0.4<£>0i Fig- 6(a) shows a broad distribution for 
( 1 1 N 1 1' )m indicating that no strict selection rules apply for 
the off-diagonal charge matrix elements in the case of fluxo- 
nium. Only at special points where the external flux in units 
of the flux quantum, Q ell i/Qo, is integer or half-integer, we 
recover an odd-even selection rule caused by the reflection 
symmetry of the potential with respect to (f = 0. Away from 
these special points, reflection symmetry is broken and strict 
selection rules disappear. Nonetheless, there is clear evidence 
that certain matrix elements are about an order of magnitude 
larger than others. The underlying regularity can be explained 
by quasi-selection rules that can be derived analytically 47 . 

In circuit QED experiments, 27 ' 29 the fluxonium circuit is 
coupled capacitively to a microwave resonator. To account 
for an additional harmonic mode observed in the experiment 48 
we include coupling to two relevant harmonic modes within a 
generalized JC Hamiltonian. The first, associated with raising 
and lowering operators a^, a, describes the fundamental mode 
of the resonator, with frequency co a . In the experiment, the 
resonator supports quarter-wavelength modes only. The low- 
est harmonic thus has a much higher frequency of 3u a , and we 
will neglect the corrections due to higher harmonics in the fol- 
lowing. The second mode represents an observed array mode 
with frequency uji, and coupling strength gt,. Both parameters 
are obtained from a fit to the spectrum in Fig. 10(a) 41 

The generic form of the effective Hamiltonian derived in the 
previous section [see Eq. (17)] now takes the concrete form 

iJ eff = w^a +w 6 b t b + ^2(ci+ki + k' 1 )\1){1\ (23) 

i 



+ r ]a ,i(ah) 2 + ?7(,;/(b t b) 2 + f^aWb 



where the energies of fluxonium levels as well as their cou- 
pling strengths are tunable by the external magnetic flux, 
= e;($ext) and g j;W = g m , ($ ext ). The ordinary ac- 
Stark shifts \a-i given in Eq. (14) are now acquire fourth- 
order contributions. Fourth-order terms are also responsible 
for the interaction terms given by ?/ a; ;(a^a) 2 | I )( 1 1 (self-Kerr) 
and £ a f, ; ;a^ab*f b| \ (cross-Kerr). These terms induce ad- 
ditional nonlinearity and result in dependence of the photon 
frequencies on the occupation numbers n a and The effec- 
tive Hamiltonian for the dispersive regime of the fluxonium 
circuit enables us to study the two central types of measure- 
ments performed in Refs. 27,29, and 30. 




FIG. 6: (Color online) (a) Magnitude of the charge matrix elements 
( I | N | I' )| between fluxonium states I and I', for external magnetic 
flux $ exl = 0.4$o- The parameters used match the experimental 
valuest. 41 Note the absence of a nearest-neighbor selection rule, (b) 
The magnitude of the charge matrix elements between ground state 
and low-lying excited states, plotted as a function of external mag- 
netic flux. An even/odd selection rule holds for zero and half-integer 
flux quantum due to parity. Away from these special flux values, sim- 
ple selection rules are absent. The vertical line marks the flux value 
0.4$ used in (a). 



A. Measurements in the dispersive regime 

The first measurement type is the direct homodyne detec- 
tion of the reflected amplitude of a single microwave tone. 
The frequency of this tone is fixed close to the bare resonator 
frequency, and the voltage of the reflected signal is recorded 
as a function of magnetic flux. This measurement primarily 
probes the dispersive shift of the resonator frequency (i.e., the 
a mode) given by \a;i to second order, and by 

fJ-a;l(n a ,nb) = E na + i : n b -i — £ ; „ ai „ i) ;/ — iO a (24) 

when including fourth order corrections. Here, the ener- 
gies E na! „, bi i are obtained as the eigenvalues of the effective 
Hamiltonian, Eq. (23). We assume that the fluxonium circuit 
occupies a fixed level I, where usually the ground state I = 
is the appropriate state maintained during the measurement. 
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Note that the Kerr terms in fourth order lead to an additional 
dependence of dispersive shifts on the excitations numbers n a 
and rib- 

Two-tone spectroscopy is the second measurement type and 
probes the fiuxonium transition frequencies via the I depen- 
dence of the resonator shift. The first tone with frequency w^i 
close to the bare resonator frequency acts in a way similar to 
the one used in the first measurement type. (As an alternative 
to the reflected amplitude, the phase shift of the reflected 
tone may be recorded.) A second drive tone is applied and 
its frequency varied over a wide range with the goal of 
inducing Rabi oscillations between the ground state and an 
excited state of the fiuxonium circuit. The transfer of proba- 
bility amplitude to higher levels is accompanied by a change 
in the dispersive shift of the resonator. In the simplest case, 
spectroscopy thus probes the corresponding change in the dis- 
persive shift given by Xa-i - Xa-o and /i a;/ - /i a;0 in second 
and fourth order, respectively. 

The detected phase shift 0; has a characteristic dependence 
on the detuning between drive frequency ojd and the effective 
resonator frequency uj a + Xa-l ( or /V')- To leading order, it 
follows the characteristic form 29 



9i 2 arctan 



2Q 



Xa-l) 



(25) 



Here, Q is the quality factor of the resonator and x ma y be 
replaced by fj, when considering fourth-order corrections. 

In our subsequent discussion of direct homodyne detec- 
tion and spectroscopy, the pole structure of the shifts Xa;i and 
fi a; l will play a crucial role. Poles are associated with res- 
onances between fiuxonium transitions and harmonic mode 
excitations, and signal the breakdown of perturbation theory 
within some frequency window centered at the pole. Since 
fiuxonium levels are tuned by varying the external magnetic 
flux $ ext , frequency windows will correspond to flux windows 
in the experiments to be discussed next. 



B. Direct homodyne detection of reflected signal 

The specific quantity observed in direct homodyne detec- 
tion is the quadrature voltage of the reflected signal. 30 This 
voltage can be expressed in terms of the phase shift via 

V = Acos(0 + O ) + C ks Acos{Bxaa + 0o ) + C. (26) 

Here, 9$ and 9' are offset phase shifts which may be present 
in experimental data, A is the amplitude of the reflected sig- 
nal and C a constant voltage offset. The approximation 
in the second step of Eq. (26) is obtained by Taylor ex- 
pansion of Eq. (25) for the drive close to resonance, i.e., 
|w<j - w tt - Xq; | < u a . 

We begin our comparison of theory with the experimental 
data from Ref. 27 with the expressions obtained in second- 
order perturbation theory. The predicted dispersive shift of 
the resonator Xa-o (while maintaining the I = ground state) 
is shown in Fig. 7(a). Its magnitude is of order 1 to 10 MHz 
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FIG. 7: (Color online) (a) Dispersive shift \a;0 of the resonator fre- 
quency for fiuxonium ground state, obtained by second-order pertur- 
bation theory. For convenient comparison with panel (b), the neg- 
ative shift — Xa-,0 is shown, (b) Comparison between experimen- 
tal data for direct homodyne detection of the reflected amplitude 
(red/light gray curve; data as published in Ref. 27) with theory fit ac- 
cording to Eq. (26) (blue/dark gray curve). The fit parameters are de- 
termined as: A = 0.0094, B = 0.0887GHz" 1 , C = 10.3052, O = 
-4.543. 



except in the immediate vicinity of the poles where the fiux- 
onium 0-1 transition crosses the resonator frequency (at flux 
values $ext/*o ~ ±0.06). 

Using Eq. (26), the dispersive shift is converted to homo- 
dyne voltage. We adjust the parameters A, B, C and 9q to 
minimize the mean-square deviations over the magnetic flux 
range ^ext/'&o € [—0.3, 0], again assuming occupation of the 
fiuxonium ground state only. We compare the resulting fit 
with the experimental data in Fig. 7(b) and find good agree- 
ment in the mentioned flux range with the exception of the 
small peak-dip structure at ^ext/'&o ~ ±0.15. Fourth-order 
corrections considered below will account for this feature. 
More significant deviations occur in the flux ranges close to 
half-integer $ e xt/$o- As we will see, fourth-order corrections 
from the effective Hamiltonian (23) alone do not lead to a sat- 
isfactory resolution, and we will discuss possible culprits for 
the persistence of deviations in this region. 

As one cause for deviations, we note that the fiuxonium 0-1 
transition reaches a minimum frequency close to 300 MHz at 
half-integer flux [see Fig. 8(d)]. For a typical temperature of 
T = 20 mK = 0.42 h GHz/fcg, thermal excitation of the I = 1 
level indeed becomes relevant. (Thermal occupation of higher 
states I > 1 remains negligibly small.) We account for thermal 
excitation under the simplifying assumption that the measure- 
ment probes a statistical mixture of the lowest two fiuxonium 
states with simple Boltzmann weights. In this case, the ther- 
mally averaged value of the dispersive shift of the resonator is 
given by 



(fJ>a;i(n ai n b )} 



Ha-o{n a ,n b ) + e eio/kBT /J, a -i(n a ,n b ) 

I _|_ g— eio/k B T 



(27) 

which we use for the remainder of this subsection. We ex- 
pect thermal effects only to be significant for the flux range 
0.35 < |3> ex t/^o| < 0.5; outside this range, eio exceeds 
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FIG. 8: (Color online) Homodyne voltage of a reflected microwave 
tone: comparison between experimental data from Ref. 27 and fourth- 
order theory for dispersive shifts, thermal averaging over fiuxonium 
levels according to Eq. (27) included (T = 20 mK). (a) Theory 
curves show dispersive shifts of the resonator frequency for occupa- 
tion numbers n a — nt — and n a = 1, n& = as indicated. The 
position of the additional pole at $ext/<&o ~ ±0.15 is in good agree- 
ment with that of the peak-dip structure observed in the experiment, 
(b) For mean photon number n a = 0.05 the theory prediction based 
on Eq. (29) also reasonably matches the amplitude of the peak-dip 
feature, (c) Additional poles are predicted if the b mode (associated 
with the array) is occupied. In each panel (a)-(c), significant devi- 
ations persist close to half-integer $ C xt/$o- (d) Fiuxonium transi- 
tion and harmonic mode frequencies, as specified by labels. Vertical 
dashed lines show alignment of resonances with corresponding poles 
in dispersive shifts and experimental features. 



2 GHz and thermal excitations of the fiuxonium device should 
be negligible. 

We now turn to the discussion of fourth-order corrections 
to the dispersive shifts. In Fig. 8(a), (b) and (c), we compare 
the same experimental data for the homodyne signal with the 
theoretical calculations now including all fourth-order correc- 
tions and taking into account thermal averaging. Specifically, 
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FIG. 9: (Color online) Level diagram illustrating a two-photon vac- 
uum Rabi splitting. In this diagram, horizontal solid lines represent 
bare fiuxonium levels which are coupled via a single photon mode. 
Here, the bare states with n a — 0, I — 4 and n a = 2, I = 
are degenerate and coupled by an effective two-photon transition. 
The resulting two-photon vacuum Rabi splitting is represented by 
the dashed lines. 



we calculate the homodyne voltage [Eq. (26)] using the same 
fit parameters as above and replace \ (/•*)■ For easy refer- 
ence of resonances, panel (d) shows the fiuxonium transition 
frequencies as well as integer combinations of the harmonic 
frequencies uj a and u>b- Relevant resonances are circled and 
labeled by roman numerals. Their alignment with the corre- 
sponding poles in (a)-(c) is indicated by vertical dashed lines. 

Figure 8(a) shows the thermally averaged homodyne volt- 
age for the case of negligible harmonic mode excitations, 
n a = 0, rift = 0, and for the case of one initial excitation 
in the resonator mode, n a = l,nt = 0. This choice of occu- 
pation numbers is motivated by the estimate of (n a ) ps 0.01 
as given in Ref. 27. In the flux region < |<& e xt/$o| ^ 0.35, 
fourth-order contributions to the dispersive shift with initial 
state n a = have no significant effect as compared to the 
second order results [Fig. 7(b)]. The dispersive shift appli- 
cable to the n a = 1 state, however, does show an additional 
pole close to |<&ext/$ol ~ 0.15 labeled by {v} in Fig. 8(a). 
As indicated in panel (d), this pole occurs due to a resonance 
between two resonator photons, 2oj a , and the 0-4 fiuxonium 
transition, €40 which we illustrate in Fig. 9. Specifically, the 
pole in /i originates from a fourth-order term in the effective 
Hamiltonian, which is associated with the V+V+V^V^ path 
of the ladder diagram [Fig. 4]: 



5o;0;i5a;;i45a;4i 3 5a;i 3 



-n a (n a -l)|0)(0|. 



f-f (w + eoiJ(2w a - e 40 )(w a + enjj ' 

(28) 

Note that this term vanishes for n a = or n a = 1. Hence, 
according to Eq. (24) it contributes to /^ a; o(n a = I, rib) 
which involves photon absorption n a = 1 —> 2 but not to 
^afi{n a = 0,n>b) where photon absorption occurs as n a = 
— > 1. This fact is easily visible in the two theory curves 
shown in Fig. 8(a). The flux position of the pole is in ex- 
cellent agreement with a similar feature in the experimen- 
tal data. We note that this resonance also manifests as the 
two-photon vacuum Rabi splitting illustrated in Fig. 9, with 
splitting size 2v / 2^ / , g a ;0l'9a;l'i/^-a;0l' (see Appendix C for 



10 



derivation). As opposed to the usual vacuum Rabi splitting, 
the two-photon splitting is proportional to g\ rather than g a . 
As a direct consequence of the absence of strict selection rules 
for fluxonium, the summation of contributions from multiple 
intermediate states /' can lead to a sizable splitting. 

To compare the amplitude of the pole feature with the ex- 
perimental reflection data, we show the weighted average of 
the dispersive shifts for n a = and n a = 1 in Fig. 8(b). 
We parametrize the respective weights Pq and Pi = 1 — Pq 
(probabilities for the two initial resonator states) in terms of 
the mean photon number n a = 1 — P . The weighted average 
is given by 

fla(n a ) =P O (Ma;i(0,0)) +(l-P )(Ma;i(l,0)). (29) 

For a mean photon number of n a = 0.05 (slightly higher than 
the value 0.01 reported in Ref. 27), we find good agreement 
between theory and experiment for the amplitude of the reso- 
nance {v}. 

In summary, for the flux region < I'&ext/'&ol S 0-3 we 
find very good agreement between theoretical prediction and 
the experimental data for the dispersive shifts, including the 
positions and amplitudes of the two resonances {v} and {vi}. 
In the flux region closer to half-integer $ ext /"I>o, however, 
agreement between experimental data and theory is weaker. 
As seen in Fig. 8(a) and (b), the poles {ii} and {iii} due to 
resonances of £40 and eio with 2uj a and cjj — u a , respectively, 
do not quantitatively match the experimental data in this flux 
region, which are dominated by a pronounced minimum at 
l^ext/^ol ~ 0.38. The resonance features in the experimen- 
tal data near |$ ext /^o| ~ 0.325 and 0.44 are absent in the 
calculated dispersive shifts of panels (a) and (b). Vice versa, 
the pole {i} predicted by theory has no correspondence in the 
experimental data. We next discuss the effects of b mode oc- 
cupations, which may give partial explanations for some of 
the mismatches. 

In panels (a) and (b) of Fig. 8, amplitudes predicted for the 
resonances {ii} and {iii} are dramatically smaller than fea- 
tures observed in the experiment at corresponding flux val- 
ues. While thermal excitation of the b mode can be ruled 
out due to the large gap between the resonant frequency 
Wf,/27r = 10.79 GHz and the frequency 0.42 GHz associated 
with a temperature of T ~ 20 mK, it is instructive to inspect 
the effects of nonequilibrium excitations of this mode. For 
this purpose, Fig. 8(c) shows the homodyne voltage in the 
presence of one b mode excitation, as obtained from the dis- 
persive shift (ii a .[(0, rib = 1)). As an important result of this 
excitation, the amplitude of resonance {iii} is amplified sig- 
nificantly. The enhancement of the resonance stems from four 
terms associated with the dual-mode ladder diagrams in Fig. 
5(b) and (c). As an example, we give the expression for one 
of them: 

gafihgb-higb-llaga-laoi^a + l)^b , qWq| 

f-f {-u> a + e i 3 )(cu b -u a — e w )(-uj a + eojj 
'1,13 

(30) 

which corresponds to the V~V b + V b V a + path. The other three 
terms arise from the paths V b + V~ V+ V b , V b + V~ V b V+ , and 
V~ V b + V a + V b . They lead to the contributions to the effective 



Hamiltonian of the form (n Q + l)n&| 0)(0 | and involve the 
denominator (u)b — u) a — eio). All of these terms indeed vanish 
for ni, — and hence do not contribute to the previous results 
in panels (a) and (b). It is thus possible that nonequilibrium 
b mode excitations are partly responsible for the deviations 
between theory and experiment in the half-integer flux region. 

By a similar mechanism, b mode excitations also lead to an 
additional predicted resonance {iv} at a flux position fairly 
close to the resonance feature in the experimental data at 
l^ext/^ol ~ 0.325. According to theory, this resonance oc- 
curs whenever the array mode excitation ujj, matches the flux- 
onium 1-3 transition e 31 . The corresponding term in the ef- 
fective Hamiltonian is 

E 2ga;13ga;3l 2 gb;l2hgb\hl / _i_1\„M\/1l «1\ 
7 . s 7 rK + W 1 M. ( 31 ) 
(-u a + eu 3 )eu 2 {uj b - e 3 i) 

«2,'3 

associated with the V b + V b V+ V~ and V+ V~ V b + V b paths in 
the dual-mode ladder diagram of Fig. 5(a). We have verified 
that a further increase in the b photon number beyond n b = 1 
[shown in Fig. 8(c)] does not improve the fit. 

In summary, we find that nonequilibrium array-mode ex- 
citations may produce significant changes in the dispersive 
shifts, some of which may point to resonance features ob- 
served in the experiment. Without a detailed understanding of 
the underlying nonequilibrium distribution and its dependence 
on magnetic flux, however, it is difficult to assess whether 
this explanation could ultimately give a quantitative match 
or whether additional array degrees of freedom, breakdown 
of perturbation theory, or dynamical effects under continuous 
driving are responsible for the observed deviations. 



C. Two-tone spectroscopy 

As explained in subsection III A, two-tone spectroscopy 
probes the change in the dispersive resonator shift due to in- 
duced Rabi oscillations between the fluxonium ground state 
and another fluxonium state (level I). When detected via the 
change in the phase of the reflected homodyne signal, the 
relevant observable is given by dio = 9i — 8q as obtained 
from Eq. (25). This phase difference is encoded by the color 
scale in Fig. 10(a), which shows experimental data from Refs. 
27,29,30. The observed transition lines correspond to the flux- 
onium 0-1, 0-2, and 0-3 transitions and to an array -mode with 
frequency cj(, = 10.79 GHz. 

Interestingly, the spectroscopy data does not only reveal the 
frequencies of relevant transitions but also contains several 
distinct phase changes along these transition lines. In par- 
ticular, the phase changes observed in the experiment can be 
classified into four types. We use the 0-2 transition line for 
illustrating these types and refer to the labels and magnified 
insets in Fig. 10(a): type I is an abrupt red-blue change, cor- 
responding to a sudden jump from positive to negative phase 
difference d^; type II is a gradual blue-white-red change, 
corresponding to a continuous change of the phase differ- 
ence from negative to positive values (or vice versa); type 
III is a gradual blue-white-blue change, corresponding to a 
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FIG. 10: (Color online) (a) Two-tone spectrum. The vertical axis represents the frequency of the sweeping tone, uja, and the horizontal axis 
represents the external magnetic flux. The color-scale represents the phase difference between the ground state and the state I excited by the 
sweeping tone, namely 6i — 6o, where red represents positive value and blue negative. The three types of color changing, I, II and III are 
illustrated by magnified insets. The color change marked by the dashed box corresponds to a spurious phase wrapping (see text), (b) Plot 
of the nonlinear dispersive shifts /-i a \i(n a , rib) — fi« ; o(na, rib), where rib = for all three curves, representing shifts with photon numbers 
n a = 0, 4, 9 respectively, (c) Plot of the nonlinear dispersive shifts Ha;2{n a ,rib) — li a ;o{n a , rib) with photon numbers n a and rib as in (b). 
(d) Plot of qudit transition frequencies and harmonic mode frequencies. Horizontal dashed lines show harmonic frequencies. Solid curves 
represent qudit frequency differences and coloring distinguishes transitions starting from 1=0, 1, and 2. Circles and squares mark zero points, 
resonances (and corresponding poles) and their alignment with experimental data. 
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the phase difference approaching zero and recovering without 
changing its sign. In principle, red-white-red color changes 

may also occur (type IV) but are not realized in this instance 
of data. 

In the following, we show that fourth-order contributions to 
the dispersive shifts explain these phase changes. The relation 
between the dispersive shift and the phase difference is given 
by 



(2Q 

9 m =2 arctan — [u dl - u a - jU a;i 



( 2Q 

2 arctan ( — [u dl - u a - fi a ;o]\- (32) 

This expression follows from Eq. (25) when setting u> d = uJdi- 
Our discussion of the different types of phase changes will 
be based on identifying magnetic flux values where the phase 
difference vanishes, 6io = (type II and III), and flux values 
where the dispersive shifts have poles such that the phase dif- 
ference may jump from tt to — tt (type I). The first condition 
is satisfied whenever fj, a; i — fi a; Q = 0, the second condition 
whenever appropriate resonances between fluxonium transi- 
tions and harmonic modes occur. We note that in the linear 
regime of the two arctan functions, the phase difference is 
simply given by 



7/0 



(33) 



Thus, in the linear regime, the two-tone spectroscopy serves 
as a direct probe of the dispersive shifts. 

For comparison with the experimental data, we calculate 
the nonlinear shifts for the fluxonium transitions 1=0— >1 and 
1=0^2, and show the corresponding differences fx a -\ — fi a -,o 
and /i a:2 — /i a; o in Fig. 10(b) and (c), respectively. In each 
case, we consider dispersive shifts for different photon num- 
bers n a = 0, 4, 9 to illustrate the dependence of the detected 
phase response on the power of the microwave probe tone. We 
use vertical dashed lines with labels to emphasize the align- 
ment of poles and zero points with the corresponding phase 
changes. According to Eq. (33), the color-encoded sign of the 
phase difference in panel (a) thus corresponds to the sign of 
— (Ha-,1 — Ma;o) m panels (b) and (c), respectively. We compare 
theoretical and experimental results for the 0-1 and 0-2 transi- 
tions in the flux regions — 0.22$ < $ ex t < and — 0.5$ < 
$ ext < where experimental data is available. 27 ' 29 ' 30 The fol- 
lowing discussion is organized according to the three observed 
types of phase changes. 



1. Type I - abrupt red-blue phase change 

The only type I phase change in Fig. 10 is marked by the 
vertical line with label {4}, and occurs due to a pole appear- 
ing in both /i a;1 — (i a .Q and /i a;2 — /J, a -Q. This pole originates 
from the fourth order contribution to /i a ;0 given in Eq. (28) 
and corresponds to a resonance between two resonator pho- 
tons and the fluxonium 0-4 transition, 2w Q = £40, see Fig. 
10(d). (This is the same resonance that also gives rise to the 
two-photon vacuum Rabi splitting illustrated in Fig. 9.) As 



noted before, this pole only occurs when n a > 1. As a result, 
the n a — curves in panels (b) and (c) do not exhibit this 
pole, and the pole becomes more pronounced as the photon 
number increases. We thus conclude that spectroscopy of the 
fluxonium transitions, in this case, also reveals information 
about photon population in the resonator. 

Two clarifications are in order. First, the absence of a type I 
phase change in the 0-2 transition due to the pole in fi a -2 — /i Q: o 
marked by the vertical line {2} falls outside our current dis- 
cussion: at this flux value, the 0-1 transition is resonant with 
w a and the discussion of the phase shift would need to include 
the hybridization of photon and fluxonium excitation. Second, 
we note that the abrupt phase change close to $ ext = — 0.4<I>o, 
marked by a dashed box in panel (a), is not associated with a 
pole and hence not of type I. Since phase changes in the ex- 
periment are defined in the interval [— tt, tt], such additional 
phase jumps may simply occur when dispersive shifts be- 
come sufficiently large so that the magnitude of 6>;o exceeds 
tt. These phase discontinuities do not involve sign changes 
in the dispersive shifts /j a; ; — fi a; o- Indeed, the occurrence 
of such a phase discontinuity close to half-integer Q ext / $0 is 
consistent with the large magnitude of the dispersive shifts 
Ma;2 — Ma:0 predicted by theory in this region. The theoreti- 
cally predicted pole corresponding to the resonance between 
£21 and u} a is slightly on the left of the dashed box. Close 
to ^ext = — 0.41$o> experimental evidence for the avoided 
crossing of the two spectral lines (the 0-2 transition and the 
0-1 transition shifted by the photon frequency uj a ) with better 
resolution can be found in Ref. 29. 



2. Type II - gradual blue-white-red change 

Gradual phase changes of type II are present in both the 0-1 
and the 0-2 transitions and are marked by the vertical line with 
label {5}. In both cases, the dispersive shift /j, a; i — /i a; o (I = 
1, 2) smoothly crosses through zero so that the phase change 
is negative for flux values < — 0.18$o> reaches zero, and then 
assumes positive values as /j, a; i — fi a - t o approaches the pole at 
position {4}. As shown in panels (b) and (c), the precise zero- 
point crossing is, in fact, photon-number dependent. We note 
that the alignment between the predicted crossings for photon 
numbers as large as n a = 9 is not perfect. Quite likely, this 
can be attributed to the breakdown of perturbation theory in 
the immediate vicinity of poles and hence especially applies 
to the zero-crossing for fj, a; 2 — (J, a ;0- As for the absence of pole 
{4} for n a = discussed above, we also expect the crossing 
{5} to disappear if no resonator photons are present. 

An additional phase change of type II is present only in the 
fluxonium 0-1 transition and is marked by line {3} and can 
be interpreted in a similar manner. Again, we find alignment 
of the zero crossing with the experimental feature only for 
rather large photon numbers, see the a n a = 9" curve in panel 
(b). An experimental study of the power dependence of this 
phase change could help shed more light on the quantitative 
comparison with theory. 
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3. Type III: gradual blue-white-blue color change 

Type III phase changes occur for both the 0-1 and 0-2 tran- 
sitions, and instances are marked by the vertical dashed lines 
labeled by { 1 } and {3} respectively. In both cases, we observe 
alignment with poles occurring in the corresponding disper- 
sive shifts. A definite prediction for type III phase changes, 
however, appears difficult based on the perturbative results. 
In general, perturbation theory will break down at the position 
of the pole and will remain unreliable in a certain flux window 
in its vicinity. As a result, predictions in this case must likely 
be based on non-perturbative methods and may possibly also 
have to take into account the dynamical aspect of the two- 
tone measurement (which are beyond the scope of this paper). 
Qualitatively, the type III phase changes are at least plausible 
given that the dispersive shifts fi a; i — /i a;0 are predominantly 
positive in the direct vicinity of both poles. 

IV. CONCLUSION 

In summary, we have presented a systematic treatment of 
fourth-order corrections to the dispersive regime of circuit 
QED. Our results, developed in Section II, are valid for a 
generic system consisting of a multi-level qudit capacitively 
coupled to one or several harmonic modes, and hence apply 
to a wide class of circuit QED systems. Our treatment, in par- 
ticular, enables the description of dispersive shifts in systems 
lacking simplifying selection rules. 

We have applied our results to the concrete case of the 
fluxonium device as realized in recent experiments. 27 ' 29 ' 30 Us- 
ing numerical diagonalization, we have obtained the relevant 
charge matrix elements which confirm the lack of selection 
rules, and have incorporated them in the perturbative treat- 
ment of the dispersive regime, including corrections up to 
fourth order. The calculated dispersive shifts allow us to com- 
pare theoretical predictions with experimental data for ho- 
modyne reflection measurements and two-tone spectroscopy 
from Refs. 27,29,30. 

The absence of selection rules is an important mechanism 
for producing sizable dispersive shifts, even if the transition 
of interest is far detuned from the resonator used for readout. 
For the fluxonium system studied in Section III, our calcu- 
lations show that dispersive shifts can indeed be as large as 
10 MHz even when the corresponding 0-1 fluxonium tran- 
sition is detuned by almost 8 GHz from the resonator. The 
lack of selection rules enables a multitude of virtual transi- 
tions to contribute to the dispersive shifts. Especially if such 
higher transition frequencies match photon resonance condi- 
tions, dispersive shifts can be surprisingly large. We also note 
that the magnitudes of matrix elements are tunable with ex- 
ternal magnetic flux, which in turn leads to the tunability of 
dispersive shifts. 

Away from half-integer & ex t/&o, we find good quantitative 
agreement for the homodyne reflection data with our theory 
prediction. This agreement also includes a resonance feature 
in the data which previously remained unexplained. The pres- 
ence of this resonance indicates a small probability of a pho- 



ton occupying the resonator, so that its amplitude may be used 
for extracting the mean photon number. Close to half-integer 
$ ext /$ > even though we find tentative agreement between 
the flux positions of several resonance features in experimen- 
tal data and theory, our calculation does not give a quantita- 
tive match. We note that the flux position coincidence of res- 
onances points to nonequilibrium array-mode excitations of 
unknown origin. 

Spectroscopy data of fluxonium samples 27 ' 29 - 30 show un- 
usual phase changes along the transition lines. We have iden- 
tified three different types of phase changes, according to 
abrupt and gradual variations of the phase with or without sign 
change. Our calculations show that these phase changes are 
closely related to poles and zero points in the dispersive shifts 
and that their occurrence may sensitively depend on photon 
numbers in the resonator. Experimental studies of the power 
dependence of spectroscopy measurements are an interesting 
subject for future study and may shed additional light on the 
origin of quantitative deviations between experiment and the- 
ory. Our results for spectroscopic phase changes show gener- 
ally good agreement for the flux positions of such resonances. 
The prediction of the specific type of the phase change re- 
mains challenging since perturbative calculations break down 
at the positions where resonances occur. Nonperturbative cal- 
culations, and taking into account the dynamics of the mea- 
surement protocol in a Master equation description, may be 
necessary to obtain such type of predictions and to improve 
quantitative agreement. An additional source of quantitative 
deviations may lie in the presence of additional array modes 
not included in our description. The spectroscopy data indeed 
shows additional levels, especially close to |$ ext /<i>o| — h< 
which warrant further investigation. 

For both types of experiments, we have identified a two- 
photon resonance near ±0.14$ ext /<I>o which manifests in the 
dispersive shifts in fourth order of perturbation theory and 
which should also lead to a two-photon vacuum Rabi split- 
ting. Experimental verification would involve tuning the drive 
frequency Ud close to the £40 transition and directly observing 
the level splitting. (In previous experiments the e 40 transition 
was outside the measured frequency range.) Further experi- 
mental verification could be achieved by detecting the corre- 
lated emission of photon pairs under vacuum Rabi oscillation. 

The accumulation of contributions to dispersive shifts in the 
absence of selection rules does not only affect ac-Stark shifts 
but can, similarly, lead to surprisingly large self-Kerr and 
cross-Kerr coefficients in fourth order which are tunable with 
magnetic flux. Making photon-photon interaction terms large 
while keeping the fundamental qudit transition off resonance, 
is particularly appealing for circuit QED lattices which have 
been discussed as quantum simulator architecture 42 ' 43 and for 
which the first experimental realizations are now underway. 44 
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Appendix A: Adiabatic elimination 

We consider a quantum mechanical system described by 
a Hilbert space H, which is composed of several subspaces, 
H = Hi © H2 © • • • For our following treatment, it is conve- 
nient to introduce the projectors 



a, J 



substituting into Eq. (A7) and consistently collecting terms up 
to A 2 , one obtains 

H' =H + [iS, H] + hiS, [iS, H}} + 0(X 3 ) 

=H + {XV + [i\Si,H Q ]} + {[i\ 2 S 2 ,H ] + [i\Si,XV] 



+ -[iXS 1 ,[iXS 1 ,H a }}} + 0(X 3 ). 



(A9) 



After imposing conditions (A4) and (A5), the resulting ex- 
pression for H' up to (and including) second order is given 
by 39 



(Al) H' = H + V di ' dS + 1^2 p a [iSi,V nd ]P a + 0(X 3 ). (A10) 



onto these subspaces. Here, the index j enumerates the eigen- 
states within each subspace H a . The projectors satisfy the 
completeness condition 



(A2) 



The goal of adiabatic elimination is to construct a Hamil- 
tonian H ' with the same eigenenergies as the original Hamil- 
tonian but without coupling between the specified subspaces. 
The elimination is facilitated by a unitary transformation, 



H' = e lb He 



iS rr-iS 



(A3) 



where the generator S is chosen such that 

P a H'P p = P a e lS He- lS P = (a ^ /3), (A4) 

which eliminates all coupling between subspaces. Condition 
(A4) only determines the off-diagonal elements of H' . To 
remove the remaining ambiguity, one imposes the constraint 
that the diagonal elements of the generator S be zero, 



P a SP a = 0. 



(A5) 



Conditions (A4) and (A5) then uniquely determine S and H' . 

To construct S and H' perturbatively and order by order, 
we decompose S into a series, 



(A6) 



where A is an auxiliary parameter for order counting. Using 
Hadamard's lemma, we rewrite 



1 1 
H' =H + [iS, H] + -[iS, [iS, H}} + ■ ■ ■ = -fiS, H] n . 

n=0 



Decomposing H' in a similar series 

oo 

H' = H + X n H' n , 



(A7) 



(A8) 



Here, the diagonal components of V are V dmg = J? P a VP a , 
and the off-diagonal elements V nd — V — T/ dlag . To first order, 
the generator is obtained as 



PaiSiPa = 



P«VP P 



E n — E 







(All) 



We then apply Eq. (A10) to a system composed of a multi- 
level qudit coupled to harmonic modes (such as the ones real- 
ized by a resonator). The non-interacting and interaction con- 
tributions to the Hamiltonian, H and V, are given by Eqs. 
(1) and (2). By our choice, the individual subspaces are one- 
dimensional each and spanned by one of the eigenstates of 
Hn, i.e., | nl)o. The integer components of n = (rij) denote 
the photon numbers of the harmonic modes (indexed by j), 
and I = 0, 1, . . . enumerates the qudit levels. 

Since V involves an increase or decrease in photon number, 
its diagonal part vanishes, V Ams = 0, and thus one finds V = 
V nd . From Eq. (All), we obtain the explicit form of the first- 
order generator 



Si 



9r,w 



9j;l'l 



,t 



en* — ujj ei'i — ujj 



l) 00 {l'\, (A12) 



where ei>i = e;/ — e; is a transition energy. Finally, Eq. (A10) 
yields the following second-order expression for H'; 

H' =flo + EE^'i"' ~ Xj;l'l)Ojaj + Xj;tl']\l)oo{l\- 
W j 

(A13) 

Here, we have abbreviated the partial dispersive shifts by 



Xj;W 



\9r,W 



\9j;W 



A 



(A14) 



By construction, H' [Eq. (A 13)] is turned into the effective 
Hamiltonian H e g in Eq. (10) by replacing bare operators by 
dressed ones, a,j — > 3j etc. The expression for the second- 
order shifts in Eq. (A14) agrees with the one derived in the 
main text [Eq. (14)]. 
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Appendix B: Numerical diagonalization of the fluxonium 
Hamiltonian 

Numerical diagonalization of the fluxonium Hamiltonian 
Hf, Eq. (22) can be achieved by using the harmonic oscilla- 
tor basis {| m)} which diagonalizes i?f|_E J= o- The oscillator 
states | m ) are generally not good approximations to the flux- 
onium eigenstates. Thus, proper care must be taken to ensure 
convergence of results with respect to the cutoff in m. Con- 
sidering ip as the position coordinate, tpo = (SEc/El) 1 / 4 
plays the role of an oscillator length. The corresponding os- 
cillator eigenenergies are given by e m — hujp^c + 1/2), 
where and c are the usual raising and lowering operators 
and uj p = ^/8ElEc / h denotes the oscillator frequency. 

In this basis, the Hamiltonian matrix has the form 



mm' c m 



Ej cos(27T$/$ )( m I cos tp\m'] 
£ , jsin(27r<l>/$o)("i I sin<^|m') 



Once the Hamiltonian matrix is diagonalized, Hf = 
J2i e/ 1 001' one obtains the charge matrix elements in the 
diagonalized basis, ( I | N 1 1' ), as follows. The charge opera- 
tor is rewritten in the oscillator basis via N = — 7= — (c — cD 
and thus one finds 



[ m I N I m ) 



Finally, one obtains 

= ^(l I m)(m\N I m!){m! \ l'), (B2) 



where ( I \ m ) and ( ml \ V ) are eigenvector amplitudes in the 
harmonic oscillator basis. 

Appendix C: Two-photon vacuum Rabi splitting 



The sine and cosine matrix elements can be expressed in terms 
of generalized Laguerre polynomials as 45 



I m I cos ip I m + 2p ) 



(Bl) 



( 2) V (m + 2p)} 
( m I sin ip I m + 2p + 1 



(-2)"^ 



2(m + 2p+ iy. 



where m,p € N. 

Alternatively, one may employ the direct evaluation of 
matrix- valued exponentials. Specifically, the cosine term may 
be re-expressed via 



cos((^ - <p m ) J^e^e-^™ + H.c. 



= 2 exp 



H.c, 



Algorithms for matrix exponentiation are standard in commer- 
cial as well as non-commercial software and libraries, e.g., 
Mathematica: MatrixExp [ ] , Matlab: expm ( ) , SciPy: 
expm ( ) and expm2 ( ) . 



Consider a multi-level qudit coupled to a single resonator 
mode a in a configuration where the bare states \h,n a )o and 
\h, n a + 2)o are degenerate. The leading-order coupling be- 
tween them is through a two-photon process involving the in- 
termediate states I n a + 1 )o, with I' = 0, 1, . . . We derive 
the effective two-photon coupling by adiabatic elimination of 
the intermediate states | l',n a + 1) and obtain an effective 
Hamiltonian H' a for the subspace spanned by the two states 
I li, n a ) and 1 12, n a + 2). By applying Eqs. (A10) and (A12) 
and using the projector P a for the relevant subspace, we find 



A 



H.c. 



(Cl) 



Here, the detuning is A a; ; 2 ;/ = 6f 2 j/ — uj a as defined previ- 
ously. Note that the effective coupling indeed involves the 
emission or absorption of two photons. Diagonalization of 
H' a yields the splitting 



2 V ( n a + 1 )("a + 1)9a\hl'9a\l'h 
; , A a; ; 2 /' 

for initial photon number n a . With n a = 0, l\ = 4 and l 2 = 0, 
we recover the expression for the two-photon vacuum Rabi 
splitting presented in the main text. 
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